Filters for High Rate Pulse Processing 
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We introduce a filter-construction metiiod for pulse processing that differs in two respects from 
that in standard optimal filtering, in which the average pulse shape and noise-power spectral density 
are combined to create a convolution filter for estimating pulse heights. First, the proposed filters 
are computed in the time domain, to avoid periodicity artifacts of the discrete Fourier transform, 
and second, orthogonality constraints are imposed on the filters, to reduce the filtering procedure's 
sensitivity to unknown baseline height and pulse tails. We analyze the proposed filters, predicting 
energy resolution under several scenarios, and apply the filters to high-rate pulse data from gamma- 
rays measured by a transition-edge-sensor microcalorimeter. 
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The extraction of physical quantities from noisy data 
streams is ubiquitous in the physical sciences. Examples 
include the determination of photon and particle ener- 
gies or incidence times in nuclear and particle physics. 
Raw data records are invariably filtered to extract the 
quantity of interest with the highest signal-to-noise, and 
extensive effort has gone into filter development. One 
important example is so-called "optimal filtering," for 
isolated pulses with amplitude proportional to photon 
energy. Filters constructed from the average pulse shape 
and the noise power spectral density are convolved with 
pulse records to estimate pulse amplitudes [l[ . This fil- 
ter is widely used in X-ray astrophysics [4I and direct 
searches for weakly interacting dark matter 

Here, we propose and demonstrate a novel method for 
optimal-filter construction for pulse processing. The pre- 
vious optimal filter is shown to be an example of a much 
larger class of filters that have new and useful properties. 
For example, optimal filters can be constructed that are 
orthogonal to exponential tails of prior pulses, prompted 
by the need to cope with high-rate pulse data. Many ap- 
plications of high-resolution photon spectroscopy require 
very large photon counts for accurate characterization of 
an absorption or emission spectrum across a broad en- 
ergy band. For example, isotopic analysis of nuclear ma- 
terials for treaty verification requires approximately 10^ 
photons in a spectrum between 60 keV and 260 keV to 
achieve uncertainty of 10^^ [4]. Low-temperature detec- 
tors can reach this goal, within limited collection periods, 
only through large arrays of elements operating at high 
photon count rates per element. Consequently, operation 
at high count rates is an active topic of research jShSj . 

The new framework departs from prior algorithms in 
two respects: (1) noise autocovariance is used in place of 
its mathematical dual, the noise power spectral density, 
to avoid the discrete Fourier transform (DFT) and en- 
able the construction to be entirely in the time domain; 



and (2) the filter optimization is subject to explicit con- 
straints beyond maximization of signal-to-noise ratio for 
isolated pulses, including for the filter length, orthog- 
onality to constants, and orthogonality to exponentials 
of one or more decay rates. (The method is related to 
constrained optimization in some other contexts. For ex- 
ample, a similar approach has recently been developed 
for designing matched filters for wavefront sensing (9|.) 
Orthogonality to exponentials can reduce or eliminate 
sensitivity to tails of prior pulses. With these additional 
constraints imposed, the filters suffer some loss of sen- 
sitivity for isolated pulses compared to filters optimized 
for that case, but they compensate by retaining resolu- 
tion with piled-up pulses and by avoiding DFT artifacts, 
including artificial periodicity. 

Processing procedure. Estimation of pulse amplitudes 
under standard filtering jfj, is optimized for isolated 
pulses. Each pulse is convolved with a filter and the 
maximum of the convolution (or a smoothed maximum 
as provided by a quadratic polynomial fit to several val- 
ues near the maximum) provides an estimate of the pulse 
amplitude. The filter, in principle, is constructed to min- 
imize the variance of this estimate, given a known pulse 
shape and known noise power spectrum. 

Continuous time model. We assume a signal / con- 
sists of a pair of pulses sitting on a baseline 

fit) = aos{t - to) + ais{t - ti) + b, 

where s is the pulse shape, ag and ai are the pulse ampli- 
tudes, to and ti are the pulse arrival times, with to < ti, 
and b is the baseline. A noisy signal consists of signal plus 
noise, m{t) — f(t) + rj{t), where the noise 77 is assumed 
to be a realization of a stationary stochastic process with 
a mean of zero and autocovariance 
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Discrete time model. The measurement apparatus 
obtains an approximation of m(iA) for i an integer, 
where A is the sample time spacing, as a convolution of 
m with a response function 



m{iA -~ t)w{t)dt, 



where w is an approximate (5-function centered at the 
origin with unit integral. We define fi, Si, and rji analo- 
gously. Our measurement model is then 



rrii^ fi + rji 



(1) 



In this approximate model, arrival times to = *oA, ti = 
iiA are assumed aligned with the samples, and known, 
to avoid interpolation issues. The pulse shape s = 
(so, . . . , s„, . . .)* is approximated by averaging many 
pulses to obtain the estimate s — (sq, . . ■ , Sn, . ■ ■)*, nor- 
malized so max s = 1 , and the noise autocovariance 
r = (ro, . . . , r„, . . .)*, given by the expectation 



(2) 



is approximated by averaging products of pulse- free sam- 
ples of the sensor output to obtain the estimate f = 
(fo,...,f„,...)*. 

Amplitude estimation. The standard procedure as- 
sumes oo = 0, computes the discrete convolution 
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of a given filter q = (go, . . . , q„_i)* with 
. . . , TO_i, mo, mi, . . . , the discrete convolution of q 
with . . . , s_i, Jo, si, . . . , where Si — for i < 0, and 
estimates oi as the ratio of their maximums 
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max; {q -k m) 



(3) 



maxi((j*s)i 

Estimate mean and variance. We seek the mean and 
variance of the amplitude estimate Si. We have 



E[{q-km)i] = aa-{q*s)t^i„+ai-{q*s),^i^+by2qj. (4) 
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We define i so that {q*s)i — max.i{q-ks)i. Under assump- 
tions of orthogonality to the prior tail and to constants. 



(q-ks) 
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(5) 



we have 
E[ai] 



E [max.i{q * m)i] 
maxj(g ★ s)i 

maxi E [(q * m)i 
maxi(q ★ s)i 



ai. 



(6) 
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FIG. 1. (color) Two scenarios, one with pile-up, are shown 
(top). From the pulse shape and noise autocovariance, a filter 
orthogonal to an exponential of tail decay, r = 3.2 ms, is 
computed (inset, separate vertical scales). Convolution of the 
filter with the signal yields peaks of essentially constant height 
(bottom) and nearly eliminates pile-up dependence. 



where the approximations are equalities under somewhat 
restrictive conditions. 

Toward a variance estimate, 

E [m^^jnii^k] = {a-oSi-io-j + aiSi-^i-j -I- b) 

X (aoSi_jo_fc + aiSi_ii_fe + b) + rj^k, 

where rj^k is the noise autocovariance of Now 

Var[ai] =E[ai2] -E[ai]^ 

E[maxi((7*m)j^] — E [maxi{q -k m)i]^ 
maxj(q ★ s)^^ 

maxi E[((7 ★ m)j^] — max^ E [{q ★ m)^]^ 



q^Rq 



maxi(g ★ s)^^ 
q*Rq dcf 



Var [ai] 



(7) 



where the variance estimate Var[ai] is defined to be the 
last expression, R is the nxn estimated covariance matrix 



with Rjk — 



f\j-k\, and s 



is the length n segment from s with q^s = maxi{q ★ s)^. 

Filter optimization. This expression for the variance 
of the amplitude estimate enables us to design filters that 
minimize the estimated variance. Var[ai] is minimized 
at a stationary point of the Lagrange function, 

Aiq,X)^q'Rq-X[q's-l,] 

where A is a Lagrange multiplier to ensure that the scale 
of q satisfies g*s = 1 = maxs. Setting the partial deriva- 
tives of A to zero and solving gives 



R- 



Var[ai] 



(8) 
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TABLE I. Dataset summary for four photon count rates. 
Sample spacing was 2.56 ^s, full stream data were kept, and 
records were later formed, including for two lengths shown. 
Discards resulted from insufficient time between triggered 
pulses, signal drops below baseline, and rises in pre-trigger 
or post-peak attributed to (untriggered) nuisance pulses. 



Dataset 


1 


2 


3 


4 


Duration (s) 


4249.50 


3096.93 


3182.13 


4583.17 


Pulses triggered 


5496 


6581 


17872 


60267 


Rate (Hz) 


1.29 


2.13 


5.62 


13.15 


Records: 10.24 ms with 5.12 ms pre-trigger 


Discards: 










pulse starts (>1) 


143 


201 


1088 


7610 


SQUID unlock 


42 


92 


585 


3648 


early peak 


22 


17 


62 


167 


pre-trigger rise 


8 


21 


93 


368 


post-peak rise 


5 


13 


65 


631 


97 keV (raw height) 


1095 


1286 


3213 


9607 


Records: 25.60 ms with 6.40 ms pre-trigger 


Discards: 










pulse starts (>1) 


243 


387 


2463 


17266 


SQUID unlock 


40 


90 


537 


2985 


early peak 


19 


17 


58 


138 


pre-trigger rise 


8 


23 


105 


606 


post-peak rise 


19 


36 


223 


1325 


97 keV (raw height) 


1067 


1228 


2938 


7565 



This solution depends on the choice, made above tacitly, 
of the length n of the convolution filter q. 

Extending beyond the prescription above, we optimize 
subject to stipulated constraints. Orthogonality to con- 
stants or exponentials of particular decay rates can be 
imposed by revising the Lagrange function. For orthog- 
onality to k vectors V — [vi ■ • - Vk] , 

A(g, A, 7) = q*Rq - A [q's - l] - q'V^, 

where 7 = (71, . . . , 7^)* are k additional Lagrange multi- 
pliers. The solution is 

q = R-W(y*R-W^ \i, Var[ai] = q*Rq, (9) 

where V — \s vi - ■ - Vk] and ei = (1, 0, . . . , 0)* is of length 
k + 1. 

Orthogonality to exponentials of a particular decay 
time constant enables filters to be less sensitive to tails 
of prior pulses. Fig. [T] illustrates the principle of these 
filters. Avoidance of the DFT, with an increase in filter 
computation cost that is very mild for filter lengths up 
to n ~ 10^, avoids false assumptions of signal and noise 
periodicity and yields nonperiodic filters. 

Experiment. Measurements were taken at NIST of 
photons from a ^^'^Gd source with a single transition- 
edge-sensor (TES) microcalorimeter at varied count 
rates (1.29, 2.13, 5.62, and 13.15 Hz), by placing the 
source at four different distances from the detector. Es- 
sentially all pulses were filtered; no attempt was made 
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FIG. 2. (color) Pulse spectrum is the absolute value of the 
discrete Fourier transform (DFT) of the average of pulses, 
near 97.431 keV line, from the highest-rate dataset. The noise 
spectrum is the average of the square of absolute value of DFT 
of pulse-free records of TES output. Records are 25.6 ms. 



Predicted energy resolution, isolated 97 keV pulse 
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FIG. 3. (color) Predicted resolution on an isolated pulse of 
four filters is shown. The filters, determined from average 
pulse shape and noise autocovariance (Fig. [J]), include the 
standard DFT-computed filter with lowest frequency bin set 
to zero 11] and proposed filters orthogonal to constants and 
zero, one, or two exponentials (ri = 6.0 ms, T2 = 1.5 ms). 



to selectively discard pulses to improve the energy res- 
olution. In extraction of pulse records from the data 
streams, pulses were lost principally due to onset within 
the prior pulse record and to occasional SQUID mode un- 
lock. Statistics for these measurements are summarized 
in Table |T1 The following analysis focuses on pulses near 
the 97.431 keV gamma-ray emission line of ^^'^Gd. 

The noise spectrum and, for comparison, the spectrum 
of the average pulse are plotted in Fig. [5] The noise 
spectrum and the DFT of the average pulse are used 
to compute the standard filter. The autocovariance and 
the average pulse (shown above in Fig. [T]) are used to 
compute the proposed filters and the predicted energy 
resolution of each. Fig. [3] shows predicted resolution ver- 
sus filter length for the proposed filters and the standard 
DFT-computed optimal filter, with the lowest frequency 
bin set to zero to reduce sensitivity to baseline drift. The 
standard filter and the proposed filter orthogonal to con- 
stants would agree, absent discretization and periodicity 
artifacts due to the DFT. This calculation is for isolated 
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FIG. 4. (color) Energy histograms near the 97.431 keV hne, 
from the 2.13 Hz dataset with 10.24 ms record length, are 
shown for the standard filter, the proposed filter orthogonal to 
constants, and the proposed filter orthogonal to constants and 
exponentials (r = 6 ms). Color denotes the pulse arrival time 
lag since the previous pulse, averaged over the histogram bin, 
and illustrates that filtering errors, concentrated in heavily 
piled-up pulses, are nearly eliminated by the filter orthogonal 
to exponentials. 



pulses; for piled-up pulses these two filters suffer bias 
problems that are significantly reduced by the filters or- 
thogonal to exponentials. The filter orthogonal to two 
exponentials, however, due to the additional constraint, 
suffers significant loss of sensitivity at short to moderate 
filter lengths, and is not considered further here. 

The performance of the other three filters is compared 
on measured pulses, and histograms near the 97.431 keV 
line are plotted for two different pulse rates in Fig. |4] 
and Fig. [5j Each histogram was fit with a Gaussian 
plus a constant to determine the energy resolution. The 
histogram bins are colored based on the pulse arrival- 
time lag from the previous pulse, averaged over the bin, 
demonstrating that errors in processing are due mainly 
to closely piled-up pulses and are significantly amelio- 
rated by the proposed filters orthogonal to exponentials. 
This effect is pronounced at the higher pulse rate, yield- 
ing much-enhanced peak height and reduced leakage for 
the filter orthogonal to exponentials. 

In Fig. m the output pulse rate, for the energy range 
97.431 ± 0.100 keV, and energy resolution are compared 
for all four input count rates and the three types of filter, 
for both short and long pulse records. At the highest rate 



FIG. 5. (color) Energy histograms as in Fig. |4l except from 
13.15 Hz dataset. At this higher rate, the errors of the first 
two filters are much more significant, as is the improvement 
offered by the third. 



and for short pulse records, the filter orthogonal to both 
constants and exponentials (r = 6 ms) offers 45 % higher 
output rate than the standard DFT-computed filter and 
40 % higher than the filter orthogonal to constants alone, 
at better energy resolution than either one. 

One important issue regarding the filters orthogonal 
to exponentials concerns their performance sensitivity 
to the choice of decay time constant r. The average 
pulse, for the TES microcalorimeter tested, was well- 
approximated over a 25.6 ms record by a linear combina- 
tion of four exponentials, with decay time constants r = 
0.018 ms, 0.144 ms, 0.963 ms, and 2.514 ms. If just the 
tail is fit, however, the constants increase considerably. 
It is evident, therefore, that no single decay rate is opti- 
mal for all arrival-time lags. Nevertheless, for the full set 
of Poisson-distributed arrival time lags, the performance 
of the filters is only mildly sensitive to the choice of time 
constant. For the highest-rate data with short records, 
over the range r = 3, . . . , 10 ms, the 97.431 ± 0.100 keV 
output pulse rate varied as 1.568 ± 0.044 Hz (mean and 
one standard deviation) and the energy resolution varied 
as 128.3 ± 3.2 eV, as compared with the t — 6 ms values 
of 1.640 Hz and 125.7 eV. 

Summary. The proposed filter construction method, 
differing from the standard procedure by being computed 
in the time domain and enabling filter optimization sub- 
ject to explicit length and orthogonality constraints, as- 
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FIG. 6. (color) The output pulse rate and energy resolution 
are compared across input count rates and three filter types 
(r = 6 ms), for both short (10.24 ms) and long (25.60 ms) 
pulse records. We note that the maximum output pulse rate 
is considerably lower than the corresponding raw pulse rate, 
because many raw pulses are due to spectral features other 
than the 97.431 keV line. 



sumes linear superposition of pulses and simple exponen- 
tial decay of pulse tails. Although these assumptions are 
satisfied rather imperfectly for the TES microcalorime- 
ter tested, the method yields notable improvement over 
standard filtering. Our tests also point to additional, 
more specialized options, such as filters optimized for a 
particular interval of arrival time lags since the previ- 
ous pulse. Such filters fit easily within this framework 
and underline that the new approach has implications 
for pulse processing in a broad range of applications. 
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